### This project replicates the analyses reported in Toshkov and Krouwel (2022) 'Beyond the U-Curve: Citizen Preferences on European Integration in Multidimensional Political Space' in European Union Politics 
### Part IV: This script prepares and analyzes the data from the 2019 EP elections in France
# Clean up VAA data and make variables -------------------------------------------------------------------------
## This part of the script shows the preparation of the datafile that is shared from the original datafile that is proprietary
# Load the data
fr.all <- read_sav('./data in/fr ep 2019/France EU19 AN.sav')

# Recode variables and filter
fr.all <- fr.all %>%
  mutate(lr = ptv_answer_ep19fr07_1242 - 5,
         proeu = ptv_answer_ep19fr01_1236 - 5,
         progr = ptv_answer_ep19fr04_1239 - 5,
         ep2019.vote.intent = recode (popup_answer_1028, 
                               '10' = 'Rassemblement National (RN)', 
                               '2' = 'La France Insoumise (FI)', 
                               '3' = 'Génération-s menée par Benoît Hamon',
                               '4' = 'Europe Ecologie Les Verts', 
                               '5' = 'Parti Socialiste (PS)', 
                               '6' = 'La République En Marche', 
                               '8' = 'Républicains (LR)', 
                               '1' = 'Parti Communiste Français (PCF)', 
                               '7' = 'UDI', 
                               '9' = 'Debout La France (DLF)', 
                               '34' = 'Other', '35' = 'No vote', '36' = NA_character_, '37' = NA_character_, .default = NA_character_),
         tk2017.vote = recode (popup_answer_1031, 
                               '4' = 'Fillon', '5' = 'Hamon', 
                               '7' = 'Le Pen', '8' = 'Macron', 
                               '9' = 'Mélenchon', '3' = 'Dupont-Aignan', 
                               '13' = NA_character_,  '14' = NA_character_,  '12' = 'No vote',
                               '0' = 'Other','1' = 'Other','2' = 'Other','6' = 'Other','10' = 'Other','11' = 'Other'),
         themes = case_when ((relevant_themes_17 == 1 ~ 'taxes_economy'),  
                             (relevant_themes_113 == 1 ~ 'money'),
                             (relevant_themes_114 == 1 ~ 'employment'), 
                             (relevant_themes_291 == 1 ~ 'globalization'),
                             (relevant_themes_292 == 1 ~ 'sovereignity'), 
                             (relevant_themes_293 == 1 ~ 'eu_policies'),
                             (relevant_themes_294 == 1 ~ 'eu_integration')),
         type = 'voter',
         
         eu.on = statement_answer_4942 - 3,
         eu.helps.protect.glob = statement_answer_4944 - 3,
         eu.weakens = statement_answer_4948 - 3,
         eu.min.wage = statement_answer_4961 - 3,
         increase.pension.age = statement_answer_4958 - 3,
         eu.eco.priority = statement_answer_4951 -3,
         nat.borders = statement_answer_4949 -3,
         fire.easy = statement_answer_4965 -3,
         eu.fp = statement_answer_4946 -3,
         redistribute = statement_answer_4962 -3,
         eu.bad = statement_answer_4945 -3,
         eu.infl.glob= statement_answer_4943 -3,
         reduce.public.spending = statement_answer_4960 -3,
         eu.workers = statement_answer_4947 -3,
         increase.wages = statement_answer_4957 -3,
         increase.state = statement_answer_4959 -3,
         state.econ.out = statement_answer_4964 -3,
         euro.out = statement_answer_4950 -3,
         deregulate = statement_answer_4963 -3
         ) %>%
  filter (country == 'France', amount_visits == 1,
          seconds_on_statements > 30*3, seconds_on_statements < 30*30) %>%
  select (lr:ncol(.))

fr.all$lr <- ifelse (fr.all$lr>=6, NA, fr.all$lr)
fr.all$proeu <- ifelse (fr.all$proeu>=6, NA, fr.all$proeu)
fr.all$progr <- ifelse (fr.all$progr>=6, NA, fr.all$progr)

# remove those who answered everything with Neutral
fr.all$zeroes <- fr.all %>% 
  select (eu.on:deregulate)  %>%
  abs () %>%
  rowSums(na.rm=T)

fr.all <- fr.all %>% filter (zeroes!=0) %>% select (-zeroes)

# create deductively derived positions on scales
fr.all <- fr.all %>%
  mutate (obj.lr = (fire.easy	- redistribute + increase.pension.age +  deregulate + reduce.public.spending - increase.wages - increase.state + state.econ.out)/8,
          obj.eu = (eu.on + eu.helps.protect.glob - eu.weakens + eu.min.wage - nat.borders + eu.fp - eu.bad + eu.infl.glob - eu.workers -euro.out )/10
          )

fr.all$id = seq.int(nrow(fr.all)) #index for merging later

# save
save (fr.all, file = './data out/fr2019_replication.RData')

# Load and subset the data (start replication here) ---------------------
load (file = './data out/fr2019_replication.RData')
# Factor analysis ----------------------------------------------------------------
fr.s <- fr.all %>% select(id, eu.on:deregulate) %>%
  filter(complete.cases(.))

# get polychoric covariance
tk.poly<-polychoric(select(fr.s, -id))

# determine number of factors
ev <- eigen(cor(select(fr.s, -id)))
ev.poly <- eigen(tk.poly$rho) #with the polychronic correlations
# number of factors

ap <- parallel(subject=nrow(fr.s),var=ncol(select(fr.s, -id)),rep=10,quantile=.95)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
plotnScree(nS)

# principal components (Table A14)
pca.1 <- principal(select(fr.s, -id), nfactors=3, rotate='none', cor='poly', scores=TRUE)
print(pca.1, digits=2, cut=.39, sort=T)

pca.out<-as.data.frame(unclass(pca.1$loadings))
pca.out <- pca.out %>%
  mutate_at(vars(starts_with('PC')), ~replace(., abs(.) < 0.4, 0)) %>%
  mutate_at(vars(starts_with('PC')), funs(round(., 2))) %>%
  arrange (-abs(PC1),-abs(PC2), -abs(PC3)) %>%
  mutate_at(vars(starts_with('PC')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
pca.out <- rbind (pca.out, c('Proportion Variance', paste0(round(prop.table(pca.1$values)[1:3]*100,0), '%')))
pca.out

pca.fr.table = htmlTable::htmlTable(pca.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.41 are not printed',
                                    align = c('l','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r'),
                                    header = c('Policy issue','PC1','PC2','PC3'),
                                    caption = "Table X. PCA results from France, 2019")
print(pca.fr.table, type = "html", file = './tables/pca_fr_table.html')

# factor analysis with polychoric correlations and oblique rotation (Table 4)
fa.1 <- fa (select(fr.s, -id), nfactors=3, rotate='Promax', cor='poly', scores=TRUE)
print(fa.1, digits=2, cut=.33, sort=T)

fa.out<-as.data.frame(unclass(fa.1$loadings))
fa.out <- fa.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR1),-abs(MR2), -abs(MR3)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa.out <- rbind (fa.out, c('Proportion Variance', paste0(round(fa.1$Vaccounted[2,]*100,0), '%')))

fa.out

fa.fr.table = htmlTable::htmlTable(fa.out, rnames = FALSE,
                                   col.columns = c('white','lightgrey'),
                                   tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                   align = c('l','r','r','r','r','r','r'),
                                   align.header = c('l','r','r','r','r','r','r'),
                                   header = colnames(fa.out),
                                   caption = "Table X. Factor Analysis results from France, 2019")
print(fa.fr.table, type = "html", file = './tables/fa_fr_table.html')

# Run the same factor model but only on the subset with party affiliation (Table A15)
fr.s2 <- fr.all %>%
  select (id, lr, eu.on:deregulate) %>%
  na.omit 

# factors with fa.poly and promax rotation
fa.2 <- fa (select(fr.s2, eu.on:deregulate), nfactors=3, rotate='Promax', cor='poly', scores=TRUE)
print(fa.2, digits=2, cut=.33, sort=T)

fa2.out<-as.data.frame(unclass(fa.2$loadings))
fa2.out <- fa2.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR1),-abs(MR2), -abs(MR3)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa2.out <- rbind (fa2.out, c('Proportion Variance', paste0(round(fa.2$Vaccounted[2,]*100,0), '%')))

fa2.fr.table = htmlTable::htmlTable(fa2.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                    align = c('l','r','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r','r'),
                                    header = colnames(fa2.out),
                                    caption = "Table SS7. Factor Analysis results from France, 2019")
print(fa2.fr.table, type = "html", file = './tables/fa2_fr_table.html')

# assign scores to cases
fr.s <- bind_cols(fr.s, data.frame(pca.1$scores), data.frame(fa.1$scores))

fr.all <- left_join (fr.all, select(fr.s, id, PC1:MR3), by='id')

# Subset-------------------------------------------------
fr.s <- fr.all %>% 
  filter (type=='voter') %>% 
  select ('lr','progr','proeu') %>% 
  filter(complete.cases(.))
# Self-placement analysis: 3D plots -------------------------------------------------
#prepare data
elevation.df = data.frame(x = fr.s$lr, y = fr.s$progr, z = fr.s$proeu)
elevation.loess = loess(z ~ x*y, data = elevation.df, degree=2, span=0.4)
xmin<- -5
xmax<- 5
ymin<- -5
ymax<- 5

elevation.fit = expand.grid(list(x = seq(xmin, xmax, 0.1), y = seq(ymin, ymax, 0.1)))
z = predict(elevation.loess, newdata = elevation.fit)
elevation.fit$Height = as.numeric(z)

#prepare matrix with x and y and rows/cols and z as the cell values
u1<-unique(elevation.fit$x)
u2<-unique(elevation.fit$y)
temp.data<-matrix(NA, nrow=length(u1), ncol=length(u2))

for (i in 1:length(u1)){
  for (j in 1:length(u2)){
    temp.data[i,j]<-elevation.fit$Height[elevation.fit$x==u1[i] & elevation.fit$y==u2[j]]
  }
}

ix <- seq(1, nrow(temp.data), length.out = 11)
iy <- seq(1, ncol(temp.data), length.out = 11)

# 3D plots in black and white (Figure 8)
png('./figures/fr 2019 eu3.png', width=900*2, height=700*2, res=96)
par(mfrow = c(1, 1), mar = c(0, 0, 0, 0))
persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Anti-Pro EU")
dev.off()


# 4 3D plots together (Figure A10 top)
png('./figures/fr 2019 eu6.png', width=900*2, height=700*2, res=96)
par(mfrow = c(2, 2), mar = c(3, 3, 3, 3))

persp3D(z = temp.data, theta = 55, phi = 20, facets = T, curtain = TRUE, 
        shade=0.8, colkey = FALSE, xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Anti-Pro EU", box=T, bty='b2',ticktype = "detailed")

persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Anti-Pro EU")

ribbon3D(z = temp.data[ix, ], along = "y", cex.axis = 0.8, 
         curtain = T, space = 0.7, theta = 25, phi = 5, shade=0.2, colkey = FALSE, xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Anti-Pro EU", 
         box=T, bty='b2', ticktype = "detailed")

hist3D(z = temp.data[ix,iy], theta = 55, phi = 20, shade=0.2, colkey = FALSE, box=T, bty='b2',ticktype = "detailed",
       xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Anti-Pro EU")

dev.off()

# now with the rayshader package (Figure A10, bottom)
temp.data2 = temp.data*300

raymat = ray_shade(temp.data2)
ambmat = ambient_shade(temp.data2)

png('./figures/raysahder_fr.png', width=1000, height=800, res=96)

temp.data2 %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(temp.data2), color="desert") %>%
  add_shadow(ray_shade(temp.data2,zscale=3,maxsearch = 300),0.5) %>%
  add_shadow(ambmat,0.5) %>%
  plot_3d(temp.data2,zscale=35,fov=0,theta=-197,zoom=0.75,phi=10, windowsize = c(1000,800),
          water=TRUE, waterdepth = 0, wateralpha = 0.5,watercolor = "lightblue",
          waterlinecolor = "white",waterlinealpha = 0.5)


render_label(temp.data2, x=1,y=1, z=1100, zscale=35,relativez=FALSE, text = "Left, Anti-Immigration",textsize = 1.7,linewidth = 2,freetype = FALSE,
             linecolor = 'darkred', textcolor = "darkred") 
render_label(temp.data2, x=101,y=1, z=1200, zscale=35,relativez=FALSE,
             linecolor = 'darkblue', textcolor = "darkblue", text = "Right, Anti-Immigration",textsize = 1.7,linewidth = 2,freetype = FALSE)
render_label(temp.data2, x=1,y=101, z=1400, zscale=35,relativez=FALSE,
             linecolor = 'darkred', textcolor = "darkred", text = "Left, Pro-Immigration",textsize = 1.7,linewidth = 2,freetype = FALSE)
render_label(temp.data2, x=101,y=101, z=1200, zscale=35, relativez=FALSE,
             linecolor = 'darkblue',textcolor = "darkblue", text = "Right, Pro-Immigration",textsize = 1.7,linewidth = 2,freetype = FALSE)

#render_label(temp.data2, x=26,y=101, z=1400, zscale=35, relativez=FALSE, dashed= TRUE,
#             linecolor = 'black', textcolor = "black", text = "Peak of EU support",textsize = 1.5,linewidth = 2,freetype = FALSE, offset=0)
render_label(temp.data2, x=85,y=15, z=300, zscale=35, relativez=FALSE, dashed= TRUE,
             linecolor = 'black', textcolor = "black", text = "Lake of Eurosceptics",textsize = 1.5, linewidth = 2,freetype = FALSE)

render_snapshot()
dev.off()


# 3D density plots --------------------------------------------------

# 3D histogram with the density mapped as color (Figure not included)
smooth <- MASS::kde2d(fr.s$lr, fr.s$progr, lims = c(-5, 5, -5, 5), n = 11)$z 

png('./figures/fr 2019 dens eu.png', width=900*2, height=700*2, res=96)
par(mfrow = c(2, 2), mar = c(3, 3, 3, 3))
color = rev(rainbow(10, start = 0/6, end = 4/6))

hist3D(z = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, colvar = smooth, col=color, colkey = FALSE,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Anti-Pro EU")

hist3D(z = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, colvar = smooth, col=color,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Anti-Pro EU")

hist3D(colvar = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, z = smooth, col=color, colkey = FALSE,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Density")

hist3D(colvar = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, z = smooth, col=color,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Anti-Pro Immigration", zlab="Density")
dev.off()

# Subjective placement and objective positions: correlation plots (Figure 7) ----------------------------------
png('./figures/fr 2019 corplot2.png', width=900, height=900, res=96)

ggpairs(
  fr.all [fr.all$type=='voter', c('obj.lr', 'obj.eu', 'lr','progr','proeu')],
  upper = list(continuous = "cor"),
  #upper = list(continuous = "density"),
  lower = list(continuous = "trends"),
  diag  = list(continuous = "barDiag"),
  axisLabels = 'show',
  columnLabels = c('Obj. Left-Right','Obj. Anti-Pro EU',
                   'Left-Right','Anti-Pro Immigration','Anti-Pro EU')) + 
  theme_bw()
dev.off()

### THE END


writeLines(capture.output(sessionInfo()), "sessionInfo.txt")
